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ABSTRACT 

Two dimensional hydrodynamical simulations of convective oxygen burning shell in the presupernova evolution 
of a 2QMq star are extended to later times. We used the VU LCAN code to simulate longer evolution times than 



previously possible. Our results confirm the previous work of Bazan & Arnett (1998) over their time span (400 s). 
However, at 1200 s, we could identify a new steady state that is significantly different than the original one 
dimensional model. There is considerable overshooting at both the top and bottom boundaries of convection zone. 
Beyond the boundaries, the convective velocity falls off exponentially, with excitation of internal modes. The 
resulting mixing greatly affect the evolution of the simulations. Connections with other works of simulation of 
convection, in which such behavior is found in a different context, are discussed. 
Subject headings: convection — methods: numerical — nucleosynthesis — supernovae 



1. INTRODUCTION 

There are many published studies of stellar convection using 
multidimensional hydrodynamical simulations, but few deal 
with convectiv e nuclear burning occurring in stellar interior (see 
Deupree 1998| for a study of core hydrogen convection). A r 



nett (1 994) (A94) and Bazan & Arnett (1994)1, B azan & Arnett 
(1998) (BA94, BA98) have studied oxygen burning shell which 
is one of the last stages in a massive star presupernova evolu- 
tion; we extend that work. 

This is an important stage in presupernova evolution because 
in this convective region several phenomena take place. In or 
near this region: (1) most of the explosive nucleosynthesis and 
production of 56 Ni and 57 Ni occurs, (2) the "mass cut" between 
collapsed and ejected matter develops, and (3) mixing of differ- 
ent layers may happen. The standard model for treating con- 
vection in one dimensional (ID) stellar evolutionary codes, the 
mixing length theory (MLT), is usually used for modeling this 
convective region as well, even though the conditions of oxy- 
gen convective burning shell are more complicated than can be 
assumed for MLT to be valid (i.e., the flow is not strongly sub- 
sonic andjhereare both energy sources and sinks in the flow). 

In Arnett (1994)| this evolutionary stage was described, and 
the first two dimensional (2D) hydrodynamical simulation of 
this problem were presented. The 145 s time interval was sim- 
ulated by these calculations, using the PROMETHEUS code, is 
much less than the d uration of this stage (10 3 to 10 5 s, see fig- 



ures 10.5 and 10.6 in Arnett (1996)1). During t his time, convec- 
tive flow was formed. |Bazan & Arnett (1994)|, Bazan & Arnett 
(1998) have modeled the evolution over a longer time interval 
of 300 to 400 s, and noticed a mixing of composition from the 
neighboring stable region near the end of the simulations. In 
addition, 

• the flow velocities were up to 10% -20% of the sound 
speed, 

• stellar structure was not altered significantly, and 

• there were strong density and composition fluctuations 
in both space and time which did not reach a statistical 



steady state. 

The main differences between this work and that of A94, 
BA94, and BA98 were (1) the use of a different hydrodynamic 
code, and (2) the simulation of a longer evolution time. The 
same initial model, equation of state, and nuclear reaction algo- 
rithms were used. 

2. THE NUMERICAL SCHEME 

For this s tudy, we use d a version of the hydrodynamical code 
VULCAN (Livne 1993). This code uses an algorithm that be- 



gins with a hydrodynamic time step, and is completed by a 
relaxation of the numerical grid, thus giving an Arbitrary La- 
grangian Eulerian (ALE) scheme. In the hydrodynamic time 
step, the Lagrangian hydrodynamic equations in two dimen- 
sions are solved explicitly or implicitly, allowing longer time 
steps to be taken if the flow is strongly subsonic. The mesh re- 
laxation phase is necessary to eliminate distortion of the cells, 
especially for flow with vorticity. We used a quasi- ID La- 
grangian relaxation, in which each ra dial row of cells kept a 
constant mass. The code was used by Glasner & Livne (1995) 



for computing convective novae outbursts, and by A sida & 
Tuchman (1997) and |Asida (2000)| to simulate convection in a 
red giant envelope. One of the adaptation made by Asida (2000) 
for simula rions of convection was the inclusion of Smagorin 
sky (1963) like sub grid scale mixing (SGSM) model, as was 
done in many two dimensional and three dimensional numeri- 
cal studies of convection. 

The equation of state and the thermonuclear reactions algo- 
rithms were essentially the same as in A94 and BA98. The 
equation of state was the sum of components for electrons, 
ions, and radiation. The nuclear reaction network used twelve 
species for helium, carbon, neon and oxygen burning. Neutrino 
cooling was included; see the above references for details. The 
initial model was the same as in A94 and BA98. The com- 
putational domain corresponds to a region containing the entire 
oxygen-burning shell in a star of 20M Q , having an initial metal- 
licity of 0.007 (about one third solar). This shell boundaries are 
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FIG. 1 . — Two dimensional velocity field for: a - 75 s, b - 150 s. The upper boundary of the oxygen shell is presented by the thick solid line in a. 



at radii of 3 x 10 8 cm and 3 x 10 9 cm that are equivalent to mass 
coordinates of 1.4M and 2.4M Q . 

We performed several simulations that were different in: the 
computational domain, spacial resolution, the initial tempera- 
ture gradient and numerical parameters of the simulation. The 
computational domain typically includes all of the oxygen shell 
as well as the neighboring layers (the inner radius was located 
at 2 x 10 8 cm in some of the simulations, and the outer radius 
was located at 8 x 10 10 cm in other simulations). The angular 
extent of the wedge ran usually from 0.35 tt to 0.65 tt radians. 

Typically, the oxygent shell was divided to ~ 120 radial 
zones, and the spacing of zones was logarithmic in radius and 
linear in angular direction (i.e., dr = r dff) with 60 angular 
zones. Rotational symmetry was assumed. These characteris- 
tics are similar to the medium resolution simulations of BA98. 
Because of small differences in the equation of state in the ini- 
tial one dimensional model and in the two dimensional code, as 
well as different radial zoning, we had two sets of simulations: 
in the first, the ID model (with interpolation) was used as it 
is, and in the second it was slightly adapted, so that the tem- 
perature gradient would be superadiabatic in all zones of the 
convection shell. 

The velocities on the inner boundary were set to be zero for 
the whole simulation, so that the inner core was a hard sphere. 
At the upper boundary there was no limitation on the veloci- 
ties, but in order to eliminate mass flow out of the computa- 
tional domain, an average radius was used to follow expansion 
or shrinking of the outer boundary. We used reflective bound- 
ary conditions on the sides (though such conditions enforce a 
downflow or upflow on the side boundaries, it was found in 
BA98 not to be important). 

3. RESULTS OF THE SIMULATIONS 

We present the results of one "standard" simulation with 172 
radial zones and 60 angular zones. The inner boundary was at 



2 x 10 8 cm and the outer was at 48 x 10 8 cm. The temperature 
gradient was slightly super adiabatic in the oxygen shell, and no 
SGSM terms were used. Most of the results of the other simu- 
lations were similar and the differences are discussed mainly in 
the end of this section. 

3.1. General Evolution of the Flow 

In our simulations, the initial velocities were zero, and the 
convective flow developed as a result of the instability from 
round off errors. Figure 1 presents the velocity field in the be- 
ginning of the simulations for times (a) 75 s and 150 s. As we 
can see, the convective flow starts at the bottom of the convec- 
tion region (i.e., the burning layer), and then moves up with 
increasing eddy size. By time 150 seconds the convective flow 
penetrates the upper boundary of the convection region (seen 
as a thick line in panel a). As a result, there is a downflow 
of carbon-rich material. This can be seen in Figure 2, which 
presents contours of carbon nucleon fraction. Panels a-e repre- 
sent times of 75, 150, 300, 600 and 1200 s. We can see that the 
downflow (panel b) penetrates the whole convective region, and 
results in mixing of carbon in this region. From comparison of 
panels c, d and e we can see that carbon abundance becomes 
more uniform, as the simulation evolves. 

This penetration of carbon is almost identical to that seen by 
BA98; see their Figure 5. The two independent hydrocodes give 
consistent results over the whole time spanned (400 seconds) by 
BA98 simulations. The small difference in the timescale for the 
carbon penetration is due to the small differences in the extent 
of the super adiabatic gradient of the initial ID model (when 
we used the initial ID model as it is, without modifying the 
temperature gradient, we got very similar results with a longer 
timescale of the penetration). 
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FIG. 2. — Carbon abundance contours for: a - 75 s, b - 150 s, c - 300 s, d - 600 s, e - 1200 s. 



3.2. Carbon Enrichment 

Carbon penetration, as well as other processes in the simula- 
tions, can be visualized by examining one dimensional averages 
of various parameters, which in principle would be equivalent 
to results from ID simulations. In Figure 3 carbon abundance 
(nucleon fraction) is plotted as a function of stellar mass coor- 
dinate. From the initial profile (solid line) we can see the lo- 
cation of the oxygen shell. At 150 s (long dashed line), carbon 
rich material has entered the convection layer, and then mixed 
through the whole region. At later times, the fluctuations in the 
abundances of carbon decrease (compare the profiles at 300 s 
-short dashed line, 600 s - dotted line and 1200 s -dotted dashed 
line), and we have a carbon nucleon fraction of ~ 5 x 10" 3 . 
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FIG. 3. — Average carbon abundance 

A change in abundance in this plot does not necessarily cor- 
responds to mixing, because this plot represents an one dimen- 
sional average of the compositions of all the material in each 
radial layer. This may be revealed by a closer look at the ap- 



parent widening of the jump in carbon abundance at the upper 
boundary of the oxygen shell at m «2.4M Q : this jump is much 
wider at 150 s than it is at later times. This "anti diffusion" is 
possible since the widening of the interface of the shell is the 
outcome of two processes: a mixture of material from the two 
sides of the interface and large scale motions that change the 
shape of the interface. Thus, at later times the interface is more 
spherical than it was at 150 s, so that averaging the abundances 
over radial layers yields a sharper jump (compare Fig 2 panels 
b and d). 
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FIG. 4. — RMS fluctuations of carbon abundance 

To better understand the carbon enrichment we present Fig- 
ure 4 which shows the RMS fluctuations of carbon nucleon 
fraction (a) as a function of stellar mass coordinate at sev- 
eral different times. The tendency toward more uniform car- 
bon mixing is clearly seen, as was indicated by Figures 3 and 
4. From this figure, we can also see, that the previously men- 
tioned widening of the oxygen shell interface is partially due 
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to changes in its shape since the fluctuations in carbon fraction 
at the boundary are relatively high, and are even higher at time 
150 s. 



Nuclear Luminosity 
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FIG. 5. — Nuclear luminosity as a function of time. 

The mixing of carbon in the burning layer causes the reaction 
rate to change significantly. In Figure 5, the nuclear luminosity 
is presented as a function of time. In the beginning there is an 
adjustment phase in which the energy production rate decreases 
to about 5 x 10 44 erg s" 1 . This transient phase is a thermal re- 
laxation due to the fact that the initial model is inconsistent: it 
does not have a two dimensional velocity field which can carry 
the convective energy flux that the ID model needed. As the 
carbon rich material penetrates to high temperature layers, it 
starts to react rapidly, yielding nuclear luminosities which in- 
crease to about hundred times the initial value. In this stage, 
the result is mainly carbon and neon consumption. Afterward, 
there is a smaller decrease in nuclear luminosity (compare time 
600 and 1200 s). This decrease is related to the small decrease 
in the average carbon abundance seen between times 600 s and 
1200 s in Figure 3. 

3.3. A Quasi-Steady State 
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FIG. 6. — Time average change in energy 

The relaxation in the carbon abundances by time 300 s to a 
slowly varying state, and the small change between times 600 
and 1200 s, suggest that the system may be close to a quasi - 
stationary state (a "steady" state on average). This one is sig- 
nificantly different from the steady state predicted by the one 
dimensional calculations used to definethe initial model. That 



the system is close to a steady state can be demonstrated by the 
following figures. Figure 6 shows the average rate of change 
in energy, as a function of mass coordinate (from time 400 to 
time 800 s). The energy produced by thermonuclear reactions 
(dashed line) is carried away by convection (solid line). The 
net change in energy is small, and balanced by changes in the 
gravitational energy (dotted line). This slow change is an ad- 
justment of the structure driven by enhanced nuclear burning 
from the carbon ingestion. There is a net heating at the bottom, 
below 1.4M Q , which will be addressed in the next section. 
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FIG. 7. — Time average change in carbon abundance 

Figure 7 shows the average rate of change in carbon nucleon 
fraction, as a function of mass coordinate (from time 400 s to 
time 800 s, and from time 800 s to time 1200 s), comparing the 
changes due to the divergence of the convective abundance flux 
(solid line) with that due to nuclear burning (long dash line). 
The nuclear consumption of carbon at the bottom is almost 
completely balanced by convective inflow; This is the burning 
zone. Over most of the convection region the change in car- 
bon abundances is small. There are larger fluctuations at the 
top (outer) boundary of the convection zone (where the carbon 
abundance increases by a factor of ~ 20, see Fig. 3). 

The decrease in nuclear luminosity (seen in Figure 5, after 
the peak at 200 s) was caused by a decrease in the net flux of 
carbon into the burning layer. This is easily seen in Figure 7 in 
the comparison of the changes in the later time interval (800 s 
to 1200 s - short dashed line) to the changes in the previous time 
interval (400 s to 800 s - solid line). The total mass of carbon 
in the oxygen convection shell is decreasing as a result of car- 
bon consumption and less mixing from the upper shells. Thus 
our convective burning is beginning to evolve on a secular time 
scale by the consumption of fuel, having approached a thermal 
"steady" state. 

in Figure 8 we present the RMS fluctuations of density (di- 
vided by the density) as a function of stellar mass coordinate 
for the standard simulation for several times. As the composi- 
tion becomes more uniform with time, the density fluctuation 
decreases throughout most of the oxygen shell to about 10~ 3 at 
time 1200 s. However at the edges of the convection zone, the 
fluctuations does not decrease and we have about 2 x 10~ 2 at 
the bottom and 7 x 10" 2 at the top, where there is a jump in 
the abundances. These values are similar to those obtained by 
BA98. 
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3.4. Interaction with Neighboring Shells 
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FIG. 9. — Temperature profile at the bottom 

One important interaction of the convective shell with neigh- 
boring shells is the mixing of carbon from the upper carbon rich 
shells, as described before. Another important interaction is 
heating of the nearest shells below the bottom of the convection 
zone: in the initial model there is a sharp decrease in temper- 
ature just below the convection zone, as convective flow starts 
to grow inside the lower parts of the convection zone there is 
some penetration of the flow to these lower temperature shells. 
As a result of this flow, the lower temperature shells are being 
heated and the higher temperature shells are being cooled. This 
heating is easily seen in Figure 9 where we present the tem- 
perature profile in the initial model, and at later times. Most 
of the heating was done by time 300 s, but even at later times, 
some heating exsisted, as can be seen in the left end of Figure 
6. Along with this "energy mixing" there is mixing of composi- 
tion in those few zones as we can see from neutron excess i] plot 
presented in Figure 10. In the initial model, there is a jump in r\ 
at the boundary of the oxygen shell (corresponds to the jump in 
composition). Due to mixing, this jump slightly moves towards 
lower mass coordinate and becomes slightly wider. From com- 
parison of Fig. 9 and 10 we see that the heating penetrates to 
deeper shells than composition changes. 



There are two important consequences of this interaction: (1) 
the heating and aditional oxygen and carbon in the few zones 
below the ID oxygen convection shell cause these zones to be 
added to the oxygen shell and, together with few zones at the 
bottom of it, to be part of the burning zone (we can see a broader 
region with temperature above 2.2 x 10 9 cm), (2) because of the 
cooling of the zones at the bottom of the oxygen shell, the tem- 
perature gradient becomes less than adiabatic in those zones 
and the convective flow exist as a result of a super adiabatic 
gradient at higher zones. 

Convective Velocity 
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FIG. 11. — Convective velocity at the bottom 

Thus we see that the neighboring stable shells greatly affect 
the oxygen convective shell. The same convective flow that 
penetrates to the neighboring shells and allow the mixing of el- 
ements and energy between the shells affects further more the 
stable shells. In Figure 1 1 we present the averaged convective 
velocity as a function of radius at the bottom of the oxygen 
shell. The dashed line corresponds to the fluctuations in radial 
component of the velocity, while the solid line corresponds to 
the amplitude of the fluctuating total velocity. The plus signs 
present points in which the pressure varies by a factor of ten. 
Three regions can be identified in this plot from right to left: 
a region of constant velocity, a region of exponential decrease, 
and another region of constant velocitis but with a noticeable 
difference between the total and the radial velocities, (i.e. the 
radial velocities are much smaller on average than the tangen- 
tial velocities) which is indicative of g modes. This g modes 
can be identified in 2D presentation of the flow (Figure 12). 
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FIG. 12. — Velocity field at the bottom. The thick line represents the com- 
position jump that indicates the boundary of the oxygen shell 

In another simulation we included a much larger shell above 
the oxygen shell. In this simulation we noticed similar flow 
characteristics of exponential decay (Figure 13) and tangential 
preference of the flow (Figure 14). From Fig. 1 1 and 13 we can 
also see that the exponential decrease at the bottom is on a scale 
which is about one third of the pressure scale height, while at 
the top it is over one pressure scale height. 



3.5. Numerical Sensitivity 

We performed several simulations with different parameters 
as mentioned in §2. As our results are consistent with those 
of BA98, and since one of the main subjects in BA98 is the 
sensitivity of the results to numerical parameters, we would not 
present here all the results from our many simulations. Instead, 
we will focus on the conclusions from those simulations: the 
most sensitive feature in the results is the amount of heating 
below the oxygen shell. When there is more heating, the tem- 
perature gradient is less than adiabatic in a larger portion of the 
(bottom of) oxygen shell. As a result, the convective flow is 
damped, carbon can not reach the burning zone and we got less 
nuclear luminosity. 

The opposite happened in simulations where we did not in- 
clude the shell below the oxygen shell: the temperature gradient 
was super adiabatic so the velocities were higher (by a factor of 
five), more carbon entered the burning zone, and the nuclear 
luminosity was higher. From comparison of simulations with 
different resolution for this case (without the bottom shell), we 
found that the results are quite robust and not sensitive to the 
numerical resolution. 

In all of our simulations the other features were essentially 
the same, namely: the evolution of the flow from the bottom, 
the penetration to neighboring shells and the enrichment of the 
oxygen shell by carbon. 
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FIG. 13. — Convective velocity at the top 
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FIG. 15. — Extra r; mixing at the bottom 




In order to check extra mixing at scales shorter than numer- 
ical resolution, we performed a simulation with SGSM mixing 
terms. This simulation started from a 2D profile of the stan- 
dard simulation at time 400 s and was carried out to time 800 s. 
The differences between this simulation and the standard sim- 
ulation were: a more uniform composition at the oxygen shell 
- the RMS fluctuations of carbon abundance at time 800 s were 
similar to the results of the standard simulation at time 1200 s, 
and significantly more mixing at the stable bottom shell. This 
can be seen in Figure 15 were we plot the neutron excess rj for 
this simulation (SGSM) and the standard simulation. In this 
plot we can see a very small difference between the two pro- 
files of the standard simulation (at time 400 s - solid line and at 
time 800 s - long dashed line). When SGSM terms were added, 
more mixing can be seen by time 600 s (short dashed line), and 
additional mixing occurs until the end of the simulation (time 
800 s - dotted line). 



FIG. 14. — Velocity field at the top. The thick line represents the compo- 
sition jump that indicates the boundary of the oxygen shell 
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4. DISCUSSION 

When we compare our results to the results of BA98 we note 
that the evolution of the simulations is similar for times reached 
previously. In both, convective flow evolves from bottom to top, 
with velocities exceeding 10% of the sound speed. The flow 
penetrates to the carbon rich region above the one dimensional 
edge of the convective region, and thus causes a significant en- 
richment of the convective region with carbon. These higher 
values of carbon abundances yield a much higher nuclear lumi- 
nosity than did the original oxygen burning. 

However, at later times, we could identify signs of a steady 
state configuration in which the average abundances changed 
very slowly, and turbulent mixing within the convective region 
was able to decrease the fluctuations in the oxygen shell. How- 
ever, density fluctuations near the boundaries of the convection 
zone (which are also composition discontinuities) did not de- 
crease at later times, and were equal to few percents. 

The most prominent feature in our simulations is the in- 
teraction of the convective shell with the neighboring shells. 
This interaction is due to penetration of the convective flow 
to the stable shells and the resulting mixing of elements and 
energy. These effects were present in all of our simulations, 
and seemed to be physically resonable and valid. However, the 
exact amount of mixing (especially at the bottom) depends on 
numerical parameters of the simulations. Moreover, since the 
initial ID model did not include such overshoot, we actually 
simulated a transient toward a more consistent model. For ex- 
ample, we mentioned that the heating of the very few neigh- 
boring zones ,at the bottom of the convection shell, lowers the 
temperature gradient below the adiabatic value. If the convec- 
tive flow is damped because of that, then the heat generation in 
these zones would evantually increase the temperature gradient, 
and convection would be restored. 

The properties of the internal modes excited in the stable 
shells are undoubtedly affected by the computational domain 
and the boundary condition. Reflective waves are probably the 
cause for the constant level of velocities we noticed at Figure 
1 1 . However, the volume below the oxygen shell in these stars 
is relatively small so a finite value of the velocity is quite plau- 
sible. 

An interesting point is the validity of the results as we per- 
formed 2D simulations (and not 3D). This is not an easy ques- 
tion to answer, however, from va rious comparisons o f 2D and 
3D simulations (see for example Kiraga et al. 1999) it seems 
that the main features of the interaction of the unstable layer 
and its neighboring layers are similar, and the overshoot and 
mixingin in both cases are comparable . As we do not claim to 
resolve this interaction quantitatively, we think that our results 
are valid. 

Penetration of the convective flow to neighboring stable re- 
gions is a common feature that exists in many multidimensional 
simulations of convection, in a varie ty of stellar problems. 
Hurlburt, Toomre, & Massaguer 1986| have studied such pen- 
etration beyond a shallow convective region in 2D, and found 
that it generated internal g modes in the stable region below the 
convection region, all the way to the bottom boundary. Two re- 
gions in this stable layer can be identified from their figure 13: 
close to the unstable layer there is a sharp decrease in the ki- 
netic energy, but it does not go to zero, rather there is an almost 
const ant level of kinetic e nergy throughout the rest of the stable 
layer. Kiraga et al. 1999 have tested penetration in both 2D and 
3D, they found that the energy fluxes in the stable layer in the 



3D are about two thirds of the 2D results and that the constant 
level is replaced with a continuas decrease when the viscos- 
ity coefficient is increased at the bottom to avoid reflection of 
waves. 

In both these studies, ideal simp lified physics was assumed. 
Freytag, Ludwig, & Steffen (1996)| studied the structure and dy- 
namics of shallow stellar surface convection zones, using two 
dimensional radiation hydrodynamic simulations with more 
"real" physics. They found convective motions extending "well 
beyond the boundary of convectively unstable region, with ver- 
tical velocities decaying exponentially with depth in the deeper 
parts of the lower overshoot region, as expected for linear g~ 
modes." They suggested to approximate the average velocity 
at the stable region as the convective velocity at the boundary 
of the unstable zone multiplied by an exponential decay factor 
with a lapse rate tha t is of order of the p ressure scale height. 

As explained by Cox (1980, §17.8) g~ modes are unstable 
within convection regions, and generally decrease exponen- 
tially with increa sing distance from the boundaries of the oscil- 
latory region. As Freytag, Ludwig, & Steffen (1996) point out, 
Landau & Lifshitz (1959) have a relevant discussion in their 
§34, where they discuss some properties of potential flow. If 
compressibility and dissipation can be neglected, a steady state 
potential flow which is periodic in some plane, must be damped 
exponentially in a direction perpendicular to that plane. Also, 
the shortest wavelengths will be damped fastest. Consequently, 
most of the overshoot is given by the largest scales which are 
driven at the interface, that is, the largest convective scale. 

We find a striking underlying unity in these results and 
ours, despite the significantly different astrophysical context in- 
volved: our simulations correspond to deep shell burning con- 
vection that is driven by nuclear burning at the bottom and 
cooled by neutrino radiation and in which spherical properties 
of the domain are important; while these studies deal with en- 
velope convection driven by photon radiation at the surface in 
plane parallel geometry. In both cases The flow extand beyond 
the formal boundaries of the convection zone, with velocities 
that decay exponentially in the stable shells producing internal 
modes. 

Near the boundaries of the convection zone, where the ve- 
locities are still large, quite efficient mixing is taking place. 
The amount of mixing beyond that region, is not entirely clear. 
Freytag, Ludwig, & Steffen (1996)| have used trajectories of test 
particles to estimate the diffusion coefficient and conclude that 
the diffusion coefficient corresponds directly to the exponen- 
tially decaying velocities. 

Some studies of stellar evolution that used these m ixing 
terms for hydrogen core convection (Herwig et al. 1997; Her- 



wig et al. 1998), found that the lapse rate of decay of the dif- 
fusion should be much smaller («2 percents) than the pressure 
scale height. They speculate that this is due to the large sta- 
bility of the layers that are near the core convection zone. Yet 
another explanation might be true, namely the flow characters 
at the stable shells and the many scales involved in mixing. 

Mixing is a complex phenomenon which take place in a very 
small length scale comparing to the length scale of the star, and 
since stars have a lower effective viscosity than in the simula- 
tions, and since they evolve for longer times, an accurate pre- 
diction of mixing is very difficult. Further, we note that motion 
does not necessarily imply mixing! In particular, oscillatory 
potential flow, which is irrotational,would involve "stretching" 
and "contracting" motions that do not give mixing, at least at 
lowest order. 
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As we can see in Figures 1 1-14, The flow in the stable shells 
is different than the flow inside the convection region, even 
though it is not purely potential flow (as revealed when examin- 
ing the vorticity). Due to these differences in the flow, it seems 
plausible that the effective mixing velocity may be less than the 
hydrodynamic velocity, and its lapse rate different as well. 

As the various shells in our model have a different compo- 
sition, mixing can be directly examined from the 2D profiles 
and the ID average of composition. However, as the scale of 
mixing is much smaller than numerical resolution, the results 
must be carefully examined. One way to examine the results is 
to compare the mixing in the standard simulation to that in the 
simulation with sub grid mixing model. There are small differ- 
ences in the average abundances in the oxygen shell between 
the two simulations, and the main difference is that with SGSM 
terms, the profile is more uniform because of the extra mixing. 
The difference between the simulations is much more promi- 
nent in the lower stable shell: the width of the abundance jump 
is much wider when SGSM terms were used, and penetrates 
to deeper zones, in which the "constant level" average velocity 
exists. 

When we used SGSM we implicitly assumed that turbulent 
flow exists from the scale of the mesh resolution to the scale 
of molecular diffusion. Since this mixing occurs in a stable re- 
gion, in which internal modes have been excited, the validity 
of this assumption is questionable. Thus it is fair to say that 
further work is needed to resolve this subject of mixing beyond 
the convective layers. 

Mazzitelli, L) Antona, & Ventura (lyyyjf have examined the 
effects of overshooting predicted by full spectrum turbulence 
model of convection (Canuto & Mazzitelli 1991, 1992) within 
the context of lithium production in AGB stars. This is an in- 
teresting example of the way stellar convection may be tested. 
A particular aspect of their work of interest here is the sugges- 
tion that symmetric overshoot may be constrained by existing 



observations. We note that a general feature of numerical sim- 
mulations is up-down asymmetry. Cooling flows are narrow, 
faster downdrafts than the corresponding updrafts, which are 
broader and slower. It is suggestive that in simulations of en- 
velope convection driven by radiative cooling at the top, there 
is an overshoot at the bottom (^Hurlburt, T bomre, & Massaguer 
1986, prey tag, Ludwig, & Steffen (1996' ), while we, who have 
simulated shell convection driven by burning at the bottom, 
have found exponential overshoot at the top (and the bottom). 
It may be that the causes of this asymmetry can be determined, 
so that a complete algorithm for the overshoot maybe devised 
for ID stellar evolutionary calculations. 

In summary, we found that stellar evolution models should 
take into account penetration and an exponential decay of the 
velocities, and the excitation of internal modes, when simu- 
lating convective burning shells. Direct hydrodynamic simula- 
tion of stellar evolution places heavy demands on computer re- 
sources, especially for carbon and neon burning stages, which 
are slower. The tendency toward steady state revealed in our 
simulations is encouraging in that simplified algorithms may 
allow the construction of an improved set of one dimensional 
models. 

We used an initial model which was derived from results 
of a one dimensional stellar evolution code, in which earlear 
stages of evolution were modeled with standard MLT Conse- 
quently, the profile we started with is unsatisfactory, and is not 
consistent with our conclusion that penetration and extra mix- 
ing should be taken into account. This should be kept in mind 
when using directly the results of our simulations. 
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